import numpy
import time
import scipy.special
import scipy.misc
import npy2cube
Параметры¶
n - основное число (1,2,3..)
l - орбитальное число (0,1,2.. n-1)
m - магнитное число (-l..+l)
Переход от декартовых координат в сферические¶
$ r (x,y,z)=\sqrt{x^2+y^2+z^2} $ - расстояние до начала координат.
$ \theta(x,y,z)=arccos\Big(\frac{z}{r(x,y,z)}\Big) $ - зенитный угол.
$ \phi(x,y,z)=arctg\Big(\frac{y}{x}\Big) $ - азимутальный угол.
Волновая функция:¶
$ a_0 $ - уменьшенный Боровский радиус.
$ R(r,n,\ell) = = \Big(\frac{2r}{na_0}\Big)^\ell e^{-\frac{r}{na_0}}L^{2\ell+1}_{n-\ell-1}\frac{2r}{na_0} $ - радиальная часть функции, где $ L_{n-\ell-1}^{2\ell+1}(...) $ - обобщённый полином Лагерра степени $ n-\ell-1 $.
$ WF(r,\theta,\phi,n,\ell,m) = R\cdot Y^m_\ell(\theta,\phi) $ - Волновая функция, заданая как произведение радиальной и угловой частей, где $ Y_{\ell}^{m}(\theta, \phi) $ - сферическая гармоника степени $ \ell $ порядка $ m $.
$ absWF(r,\theta,\phi,n,\ell,m) = |WF(n,\ell,m,r,\theta,\phi)|^2 $ - плотность вероятности нахождения частицы в данной точке пространства в данный момент времени.
def w(n,l,m,d):
x,y,z = numpy.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]
r = lambda x,y,z: numpy.sqrt(x**2+y**2+z**2)
theta = lambda x,y,z: numpy.arccos(z/r(x,y,z))
phi = lambda x,y,z: numpy.arctan(y/x)
a0 = 1.
R = lambda r,n,l: (2*r/n/a0)**l * numpy.exp(-r/n/a0) * scipy.special.genlaguerre(n-l-1,2*l+1)(2*r/n/a0)
WF = lambda r,theta,phi,n,l,m: R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
absWF = lambda r,theta,phi,n,l,m: numpy.absolute(WF(r,theta,phi,n,l,m))**2
return WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
Расчет значений для первых трех уровней. Создание Cube-файлов.¶
Также пришлось немного изменить скрипт npy2cude для совместимости python3
d = 30
step = float(2.*d/29)
# Зададим цикл по перебору квантовых чисел
for n in range(0,4):
for l in range(0,n):
for m in range(0,l+1,1):
grid = w(n,l,m,d)
name ='%s-%s-%s' % (n,l,m)
npy2cube.npy2cube(grid,(-d,-d,-d),(step,step,step),name+'.cube')
PyMOL¶
from IPython.display import display,Image
from xmlrpc.client import ServerProxy
cmd = ServerProxy(uri='http://localhost:9124/RPC2')
# Откуда эти цифры?
# vol obj -> coloring -> panel -> *manual change* -> get colors as script
# cmd.volume_ramp_new('name of profile', [\
# pointi: thickness, R, G, B, opacity,\ ...])
cmd.reinitialize()
cmd.volume_ramp_new('ramp007', [\
0.002, 0.00, 0.00, 1.00, 0.00, \
0.01, 0.00, 1.00, 1.00, 0.20, \
0.015, 0.00, 0.00, 1.00, 0.00, \
])
cmd.set('ray_volume')
cmd.load('/Users/alex/Desktop/ws/hse/golovin/task4/1-0-0.cube')
cmd.volume('vol', '1-0-0')
cmd.volume_color ('vol','ramp007')
cmd.ray(300, 300)
cmd.png('/Users/alex/Desktop/ws/hse/golovin/task4/1-0-0.png')
time.sleep(1)
Image('/Users/alex/Desktop/ws/hse/golovin/task4/1-0-0.png')
cmd.reinitialize()
time.sleep(2)
cmd.set('ray_volume')
cmd.volume_ramp_new('ramp007', [\
-0.05, 0.00, 0.00, 1.00, 0.00, \
-0.01, 0.00, 1.00, 1.00, 0.20, \
-0.015, 0.00, 0.00, 1.00, 0.00, \
0.005, 0.00, 0.00, 1.00, 0.00, \
0.01, 0.00, 1.00, 1.00, 0.20, \
0.01, 0.00, 0.00, 1.00, 0.00, \
])
names = []
for n in range(0,4):
for l in range(0,n):
for m in range(0,l+1,1):
name='%s-%s-%s' % (n,l,m)
file = '/Users/alex/Desktop/ws/hse/golovin/task4/' + name + '.cube'
cmd.load(file)
cmd.volume('vol', name)
cmd.volume_color ('vol','ramp007')
time.sleep(10)
cmd.png('/Users/alex/Desktop/ws/hse/golovin/task4/' + name + '.png')
time.sleep(5)
names.append(name)
for name in names:
display(Image('/Users/alex/Desktop/ws/hse/golovin/task4/' + name + '.png'))
ORCA¶
определенная схожесть имеется
cmd.set('ray_volume')
cmd.volume_ramp_new('ramp007', [\
-0.05, 0.00, 0.00, 1.00, 0.00, \
-0.01, 0.00, 1.00, 1.00, 0.20, \
-0.015, 0.00, 0.00, 1.00, 0.00, \
0.005, 0.00, 0.00, 1.00, 0.00, \
0.01, 0.00, 1.00, 1.00, 0.20, \
0.015, 0.00, 0.00, 1.00, 0.00, \
])
for i in range(1,5):
name='H-%s'% i
file = '/Users/alex/Desktop/ws/hse/golovin/task4/' + name + '.cube'
cmd.load(file)
cmd.volume('vol', name)
time.sleep(10)
cmd.png('/Users/alex/Desktop/ws/hse/golovin/task4/' + name + '.png')
time.sleep(5)
for n in range(1,5):
name = 'H-%s' % n
display(Image('/Users/alex/Desktop/ws/hse/golovin/task4/' + name + '.png'), name)
